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The Kotliar-Ruckenstein slave-boson scheme is used to allow for an unrestricted variation of 
the bosonic and fermionic fields on the saddle-point level. Various inhomogeneous solutions, such 
as spin polarons and domain walls, are discussed within the two-dimensional Hubbard model and 
compared with results of unrestricted Hartree-Fock (HF) calculations. We find that the present 
approach drastically reduces the polarization of these states and leads to increased delocalized wave 
functions as compared to the HF model. The interaction between two spin-polarons turns out to 
be attractive over a wide range of the on-site repulsion U. In addition we obtain the crossover from 
vertical to diagonal domain walls at a higher value of U than predicted by HF. 
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I. INTRODUCTION 



The unrestricted Hartree-Fock (HF) approach has 
turned out to be a powerful tool in the calculation of in- 
homogeneous states in the Hubbard model ■ Among 
them magnetic polarons , domain walls |g-g,0j29| 

and vortex solutions || have been extensively studied 
within the context of high T c superconductors (HTSC). 

Schrieffer, Wen and Zhang || have proposed the so- 
called spin-bag mechanism according to which a hole cou- 
ples to the spin density of the antiferromagnetically (AF) 
ordered background of the CuC>2 planes thus creating a 
local reduction of the AF order parameter. Within this 
model two holes are attracted by sharing a common bag. 
In the simplest model the resulting superconducting or- 
der parameter is of the order of the spin-density wave 
gap which is large enough to lead to high-temperature 
superconductivity. The concept of spin polaron forma- 
tion is also successful in explaining the phase diagram in 
the low and intermediate doping regime of the HTSC's. 
According to the concept of microscopic electronic phase 
separation |p|,^0|, the doping induced spin polarons (or 
spin clusters) form a conducting subsystem in an back- 
ground dominated by strong antifcrromagnetic correla- 
tions. Due to the competition between short-range at- 
traction and the long-range Coulomb repulsion, these 
polarons do not form a homogeneous sea in an antiferro- 
magnetic background. Instead, the holes are confined to 
a percolative network which dimension D is lower than 
two (1 < D < 2). The resulting inhomogeneity in the 
electronic subsystem has also important implications on 
the pairing mechansim ]l2[ ]. 

Besides spin-polarons, domain wall solutions presently 
also attract a lot of interest in the field of HTSC. 
Stripe correlations have been observed in L^-zSr^CuCU 



]l3| fl5| , in nickel oxide analogues of the copper oxides |l6| ] 
and in Lai.6-xNdo.4Sr x Cu04 (l^jl^]. In the two latter 
systems the density fluctuations are pinned by the under- 
lying lattice structure giving rise to a static charge- and 
spin-density wave. The low-temperature orthorhombic 
phase in the nickel oxides stabilizes the stripe structure 
along the diagonals whereas the low-temperature tetrago- 
nal structure in Lai.6-2;Ndo.4Sr x Cu04 favors the pinning 
of a horizontal stripe phase. Also in the Cu02-planes of 
Bi2212 compounds the existence of charge stripe order 
has been demonstrated with extended x-ray-absorption 
fine structure (EXAFS) experiments ||l9|| . 

The stripe instability was predicted theoretically in Q| 
within a HF formalism applied to the extended Hubbard 
model. Also most of the further investigations on the 
striped phases [||J|J^,^] have been carried out within 
standard HF theory. However, the solutions obtained 
with HF are poor variational wave functions since they 
are much too high in energy. Within HF, the only mech- 
anism avoiding double occupancy in order to reduce the 
Hubbard repulsion is to renormalize the spin-dependent 
on-site energy e,^ = U(rii_- a ). Therefore the commen- 
surate antiferromagnetic phase displays as an alternating 
shift of the spin-up and spin-down one-particle levels re- 
spectively, which overestimates by far the polarization of 
the AF order. 

It is well known that a part of the correlations between 
electrons of opposite spins can be accounted for by us- 
ing the Gutzwiller projector pl| , p2"| |. For the AF ordered 
system this approach leads to wave functions that are 
very close in energy to the solution obtained by Quan- 
tum Monte Carlo simulations [ p3|j2^ ]. In addition, the 
Gutzwiller wave function leads to a significantly lower 
spin polarization in the intermediate U regime as com- 
pared to the HF result p3| . The Gutzwiller approach has 
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also been used to improve the solutions of inhomogeneous 
states such as spin polarons (f^Hl and domain walls 
p5|p6| . However, due to computional limitations these 
calculations where done with an Ansatz for the charge- 
and spin-density profile, i.e the number of variational pa- 
rameters for charge and spin at the different lattice sites 
has been strongly reduced. 

In this paper we overcome the limitations discussed 
above and present results for unrestriced spin-polaronic 
and kink type Gutzwiller wave functions. Our approach 
is based on the representation of the Hubbard model in 
terms of fermions and slave bosons due to Kotliar and 
Ruckenstein (KR) (27J. The KR formulation is a func- 
tional integral method that reproduces the Gutzwiller 
solution at the saddle-point level. The advantage of this 
method is that it provides a systematic way to improve 
the solution by expanding the fields around the saddle 
point. The KR approach is a good starting point for our 
purposes, since at the saddle point level we can immedi- 
ately identify the variational parameters with the bosonic 
and fermionic fields. Since we don't make any assump- 
tion on the spatial symmetry of these fields we therefore 
obtain various inhomogeneous textures which stability 
depends on doping and the value of the Hubbard repul- 
sion U. In the present paper we will concentrate on the 
stability and shape of spin-polaronic and domain wall so- 
lutions, respectively. It turns out that for commensurate 
filling (i.e. one hole along the wall) the parameter range 
for the occurence of stripes is significantly enlarged in 
comparison to the HF calculation. Within the investi- 
gated range of U (4t < U < lOt) we don't observe a 
crossover to a polaronic Wigner crystal for commensu- 
rate doping, which means that the interaction between 
the holes keeps attractive up to very large values of U. 



The rest of the paper is organized as follows: In Sec. 

II we give a detailed description of the formalism, in Sec. 

III we present the results for spin polaronic and domain 
wall solutions respectively, and in Sec. IV we summarize 
our conclusions. 



II. MODEL AND FORMALISM 

We consider the two-dimensional Hubbard model on a 
square lattice, with hopping restricted to nearest neigh- 
bors (indicated by the bracket < i, j >) 



H 



where destroys (creates) an electron with spin a at 

site i, and n.; j(T = cj a Ci t „. U is the on-site Hubbard repul- 
sion and t the transfer parameter. For the calculations 
in Sec. Ill we take t=f. Following KR we enlarge the 
original Hilbert space by introducing four subsidiary bo- 



-t 



E 

<ij>,a 



U J2 



(1) 



son fields e 



(t) D (t) Jt) 



and (V for each site i. These 
operators stand for the annihilation (creation) of empty, 
singly occupied states with spin up or down, and dou- 
bly occupied sites, respectively. Since there are only four 
possible states per site, these boson projection operators 
must satisfy the completeness condition 



Furthermore 



dldi 



d]di 



Then, in the physical subspace defined by Eqs. (|]||) the Hamiltonian ([!]) takes the form 

H = -t z tA,* c i,* z i,° + U Y1 4 d i 



<ij>,a 



1 



1 



(2) 

(3) 

(4) 
(5) 



and has the same matrix elements as those calculated for (|j) in the original Hilbert space. 



In the saddle-point approximation, all bosonic opera- 
tors are treated as numbers. The resulting effective one- 
particle Hamiltonian describes the dynamics of particles 
with modulated hopping amplitude and can be diagonal- 
ized by the transformation 



ok 



(6) 



where the orthogonality of the transformation requires 

5>k(*)av(9) = V (7) 



Given a system with N e i particles we finally obtain for 
the total energy 



E tot = -t <^E $ U fc )^( fc ) + uJ2 d * 



k=l 



(8) 



which has to be evaluated within the constraints (^|Jj, 
^) . This is achieved by adding these constraints quadrat- 
ically to Eq. (g) 
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i!fci = Ai£(e? + J>? iir + d?-l) a ( 9a ) 
E C 2 = A 2 **iAk)$iAk) - s\a ~ dj) 2 (9b) 

i,cr k 

E C3 = A 3 ^(E KAk)^A<l) ~ 5 kq f (9c) 

k,q i,(T 

E C 4 = A 4 (^$* T (fc)$^ T (fc) - 7V T ) 2 (9d) 

k . i 

E C5 = A 5 (^ $* a (fc)$ a (fe) - A^) 2 (9e) 

We have added the last two conditions that turn out to 
be very convenient because they allow to define the total 
number of spin up and down particles Af + N | = N e j. 
The energy functional E{^i <a (k), e^, Si i(T , di} = E to t + 
Eci + Eci + Ec3 + Ec4 + Ecs now has to be minimized 
with respect to the fermionic and bosonic fields. Since 
the KR theory does not preserve the spin-rotation in- 
variance of the original hamiltonian all variational pa- 
rameters can be taken as real numbers [p8[ . For the 
minimization procedure we have used a standard conju- 
gate gradient algorithm. The gradients of the functional 
E{§i ir7 (k), ei, Si i<T , di} can be calculated analytically and 
convergence is checked by evaluating the norm of the gra- 
dient. The accuracy of the solution can be controlled by 
calculating the value of Eci + Ec2 + Ec3 + Eca + Ec5 at 
the end of the iteration procedure. We generally have set 
the values of the Lagrange parameters A1...A5 to 10 4 — 10 5 
which leads to an estimated Error at ss 0.0002. 

In principle one could start the calculation with a ran- 
dom configuration of the fermionic and bosonic fields. 
However, for a doped system there exist different self- 
consistent solutions which are close in energy and deter- 
mining the most favorable can be difficult. Therefore, we 
have generally started from the unrestricted HF solutions 
for spin polaronic or domain wall phases. The order of 
magnitude of the time needed to get convergence is half 
an hour on a SGI Indy workstation. 

III. RESULTS 

In a first step we consider a single spin-polaron and 
compare the spin and charge profiles obtained within the 
HF and unrestricted SB approximation, respectively. We 
then extend the calculations to domain wall type solu- 
tions and study the stability of these textures as a func- 
tion of U. Before this, we evaluate the ground state en- 
ergy and AF polarization for the half filled system and 



compare our calculations with the results of Yokoyama 
and Shiba p3|| . Since for this case their AF Gutzwiller 
variational approach (AFGF) is equivalent to the saddle- 
point approximation of our slave boson method we find 
perfect agreement within the numerical error. In con- 
trast to the HF scheme, the AFGF leads to a signicfi- 
cant reduction of the magnetization in the intermediate 
U regime. On the other hand, both methods converge 
rapidly to fully polarized S z components in the large U- 
limit. One should note that Quantum Monte Carlo sim- 
ulations |Q still lead to much smaller values in the large 
U limit. 

A. Spin Polarons 

The formation of spin polarons in the 2D Hubbard 
model results from the competition between kinetic en- 
ergy gain and magnetic energy loss when doping the sys- 
tem away from half filling. Let us first consider the case 
where a particle with spin down has been removed from 
the half filled, antiferromagnetically (AF) ordered lat- 
tice. If the vacancy is immobile, the cost in the magnetic 
energy is 4J^ t 2 /U in the large U limit. However, this 
vacancy can gain kinetic energy via virtual hopping pro- 
cesses to the nearest neighbor sites, thus mixing some 
probability of spin up occupation to the site where the 
particle has been removed from. Therefore, the result- 
ing spin-density profile therefore has an inverted order 
parameter at the site where the charge is located. Ac- 
cording to Nagaoka's theorem |2(]] the size of this ferro- 
magnetic core is expected to be very large for large values 
of U. 

Fig. [I] show the charge and spin density profile of such 
a polaron obtained by unrestricted HF calculations and 
with the present method, respectively. The calculation 
was done on a 8 x 8 lattice with on-site repulsion U=6t. 
In both methods, the doped hole is mainly localized at 
site (4,5), But, whereas the charge at this site is reduced 
to (n) — 0.53 within the HF approach, the unrestricted 
SB method gives a value of (n) = 0.73 only. Moreover, 
the AF order parameter Af = (— l) l ' r+ ' ly at the po- 
laron center is much less affected in the SB mean field 
(MF) treatment (Af = —0.07) than in the case of un- 
restricted HF (Af = —0.21). This discrepancy can be 
understood as follows. The HF theory only renormalizes 
the spin dependent on-site energies. The removal of a 
spin down particle at site i leads to a relaxation of the 
spin up on-site level at this site. As a consequence the 
alternating on-site level shift, describing the AF order, 
is changed at site i where 5 neighbored spin up states 
have now acquired nearly the same energy. Thus there 
is a strong hybridization between the spin-up states and 
one obtaines a large value for the reversed spin order 
parameter at the central site i. By contrast the kinetic 
energy in the spin-down channel between site i and its 
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nearest neighbors is very much reduced since the cor- 
responding spin-down on-site level is pushed to a high 
energy. Note that the HF method always leads to a very 
large spin-polarization, since this is the only way within 
this approximation to minimize the on-site repulsion. 

Let us now consider the removal of a spin down par- 
ticle at site i when calculated within the SB aproxima- 
tion. Then the spin-up hopping channel allows for the 
hybridization of neighboring spin up states with site i, 
which is comparable to the HF approach. However, since 
we now have an additional variational parameter per site 
(i.e. the boson field {di}) double occupancy can be min- 
imized at site i, whithout very much reducing the kinetic 
energy in the spin-down channel between site i and its 
nearest neighbors. In fact the reason why the Gutziller 
wave function leads to a lower energy than HF theory 
is that there one can minimize double occupancy while 
keeping the kinetic energy at a higher value as compared 
to the HF approach. In other words, the localization be- 
havior is much more pronounced in the HF treatment 
than within the SB formalism. 

The energy difference of the solutions in Fig. [I] (U=6t) 
is of the same order of magnitude than for the homoge- 
neous AF solutions of Ref. p3| . The total energy per 
site calculated within the HF approximation is E HF = 
-0.607t and for the SB method we get E SB = -0.645i. 
Fig. [2] shows the order parameter Af as a function of 
U at the perturbed sites and far from the spin polaron, 
evaluated with HF and SB approximations respectively. 
It turns out that the SB method drastically reduces the 
polarization at the center of the polaron as compared to 
HF. As can be further seen on Fig. |^, this reduction of 
polarization is much stronger than in the residual AF 
ordered plane. Within the SB scheme the spin polaron 
acquires a ferromagnetic core for values U > 5t whereas 
this limit in HF is already achieved for U ~ 3.3i. 

Finally we evaluate the static interaction between two 
spin-polarons which is of special importance with regard 
to the spin-bag model |^| and phase separation scenar- 
ios || of the high-T c superconductors. To calculate the 
binding energy we consider two holes with opposite spins 
placed at two neighbored sites and compare the energy of 
this configuration with the energy of two infinitely sep- 
arated spin-polarons. It should be noted that a positive 
binding energy means attraction. The results are plotted 
in Fig. |^. Within the HF approach we obtain an interac- 
tion between the two polarons which is attractive up to 
U rs 6.5t, in agreement with the results of Ref. Q]. The 
binding energy displays a maximum at U ~ 4.5i. How- 
ever, within the SB approximation the parameter space 
of attraction is considerably enlarged. The maximum 
of the binding energy, which is approximately twice the 
value of the HF calculation, now occurs at U « 8.5t. 
Unfortunately, convergence of our variational approach 
becomes very slow in the very large U-regime. In this 
limit one has to increase very much the parameters Xi 



in eqs. (|9[) in order to keep the constraint induced error 
within the desired limits. This fact causes the difficulty 
in exploring the whole parameter range of attraction in 
Fig. I 



B. Domain wall solutions 

For intermediate values of the on-site Hubbard repul- 
sion U Hartree-Fock theory predicts the existence of do- 
main wall solutions, where the doped charged carriers are 
localized within a stripe in horizontal or diagonal direc- 
tion. This stripe separates two AF ordered regions with 
opposite sign in the AF order parameter. Within the 
HF approximation it was shown Q that there is a tran- 
sition from horizontal to diagonal stripes when the ratio 
between on-site repulsion U and transfer integral t ex- 
ceeds the critical value of U/t ~ 3.6. From a constrained 
Gutzwiller variation of hyperbolic-type domain walls it 
was concluded in pq ], that this limit probably is shifted 
to much lower values. 

In the following we will compare energies for diagonal 
and horizontal stripes with the energy of isolated spin- 
polarons using the unrestricted SB scheme. The calcula- 
tions are made for different lattice sizes (typically 17x4, 
13 x 6, 9 x 8) by applying appropriate boundary condi- 
tions for each domain wall type (see Fig. |J). The choice 
of the supercell is, of course, a delicate issue since peri- 
odic boundary conditions in principle require an 'uneven 
x uneven' lattice for a diagonal wall and an 'uneven x 
even' lattice for a vertical stripe. To avoid the comparison 
between different lattice sizes for different domain wall 
types we have choosen the 'shifted boundaries' shown in 
Fig. ^b. This means that along the y-direction the su- 
percells are shifted by the extend of the diagonal domain 
wall in x-direction. Consider for example a diagonal do- 
main wall on a 13 x 6 lattice. Then the site (N X) 1) is 
connected with site (N x + 6,6) for N x < 8 and with site 
(N x +6-13,6) for N x > 8. Calculating the energy per site 
for the half-filled AF ordered system, we find that the two 
kinds of boundary conditions in Fig. || differ in the result 
by 0.1% for U~ 3t. This difference rapidly vanishes with 
increasing U. 

In Fig. ^ we show the charge- and spin density profile 
of a vertical domain wall calculated with unrestricted HF 
and SB approximations respectively. The doping corre- 
sponds to one hole per site in the domain wall. As for 
the spin polarons studied in the previous section the SB 
result displays a charge- and spin profile that is consid- 
erably enlarged compared to the HF solution. This is in 
agreement with the calculations of Ref. |2(J . The dip in 
the HF charge profile is nearly twice the value of the SB 
approximation. 

Fig. U shows the energy per hole as a function of U for 
various textures such as domain walls and spin polarons. 
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We also have investigated the half-filled wall with the on- 
wall quadrupling of the period which has been intensively 
discussed by Zaanen and Oles in Ref. Jp9[ , The energy 
has been calculated in a standard way |2S(| by comparing 
the energy of each texture with the energy of the ho- 
mogeneous AF ordered state with the same number of 
holes (compared to half filling). In case of the diagonal 
wall we additionally have choosen the same shift of the 
boundaries for the reference AF lattice. 

Within the HF approximation we obtain a crossover 
from vertical to diagonal domain walls for Uw 3.8i and 
a crossover to isolated spin-polarons at U~ 8t. The en- 
ergy for half-filled walls is always higher than for iso- 
lated polarons. These results are in complete agreement 
with earlier HF studies of the two-dimensional Hubbard 
model H^]. However, the range of stability for the ver- 
tical stripe solution is considerably enlarged in the unre- 
stricted SB approximation where we obtain the crossover 
at U« 5.7t This result is supported by Lancos diago- 
nalization studies of the tJ-model |5(J and Monte Carlo 
methods of the one-band Hubbard model []3l|~|33f . These 
works report a shift of the static spin structure factor 
peak in the vertical direction when the charge density 
is reduced away from half-filling in agreement with our 
findings. Also recent studies of the 2D tJ model within a 
density matrix renormalization group approach |34| are 
in agreement with a vertical striped phase. In addition 
we don't observe a crossover to a spin-polaronic Wigner 
crystal for the considered range of U (U < lOt). In- 
stead the energy of half-filled walls turns out to be lower 
than the energy of spin polarons. In fact, since in the 
SB approach the range of attraction between the holes is 
considerably enhanced, these textures are energetically 
favored which have the holes placed nearby each other. 
However, for very large values of U we again expect a 
decay of the stripe into isolated polarons. 

IV. CONCLUSIONS 

We have shown that the unrestricted SB saddle-point 
approximation is a simple and powerful tool to improve 
inhomogeneous solutions obtained by the HF method. It 
turns out that this approach leads to a strong reduction 
of the spin-polarization of these inhomogeneities in the 
intermediate U regime. The most relevant feature of the 
SB approximation, however, is the considerably enlarged 
range of attraction between spin-polarons in comparison 
with the HF method. This result has also a strong im- 
pact on the spin-bag model of high T c superconductivity, 
since HF theory has restricted the validity of this model 
to small values of the Hubbard repulsion U. 

Regarding the domain- wall phases, we find that the 
crossover from vertical to horizontal stripes is shifted to 
higher values of the on-site repulsion U than predicted 
by HF theory. Moreover we don't observe a crossover 



to isolated polarons for U< lOt. This result supports 
the description of domain wall structures in the La2Ni04 
compounds |35| where it is generally argued that the oc- 
curence of diagonal walls is a result of mean field theory. 
Within the present approach the range of stability of di- 
agonal walls is considerably enhanced in comparison to 
the HF approximation, where for U« 8t the walls decay 
into isolated polarons. 
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FIG. 1. Charge- ((n»)) and Spin- (Af ) density profiles for a spin polaron on a 8 x 8 lattice, a) HF approximation; b) SB 
approximation. The Hubbard on-site repulsion is U = 6t. 

FIG. 2. Spin-density order parameter Af at the center of the spin polaron (two lower curves) and at maximum distance 
from the polaron (two upper curves). Solid line: SB approximation; Dashed line: HF approximation. 

FIG. 3. Binding energy of a pair of polarons with opposite spins placed at neighbored sites on a 8 x 8 lattice. Dashed line: 
HF approximation; Solid line: SB approximation 

FIG. 4. Sketch showing the boundary conditions in y-direction which have been applied to describe vertical (a) and diagonal 
(b) domain walls respectively. In x-direction we have taken periodic conditions. The line crossings correspond to the lattice 
sites and the dashed lattices mark the choosen periodicity in the y-direction. 

FIG. 5. Charge- (a) and Spin- (b) density profiles in x-direction for a vertical domain wall on a 17 x 4 lattice. The number of 
holes is 4 since the best energy is obtained when there is one hole per site in the wall. Solid lines: SB approximation; Dashed 
lines: HF approximation; Short dashed: cosh- (for charge) and tanh- (for spin) functional fit to the SB solution. 



FIG. 6. Binding energy per hole for vertical stripes (full lines), and diagonal stripes (dashed lines) with one hole per site 
along the wall, vertical half-filled walls (dotted lines) and isolated spin polarons (dashed-dotted lines). The binding energy is 
defined as the difference in energy between a given texture and the homogeneous AF ordered lattice with the same doping. The 
domain wall solutions have been calculated on a 9 x 8 lattice, polarons on a 8 x 8 lattice. a)HF approach; b)SB approximation. 
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